rm(list=ls()) #clear

library(psych)
library(ggplot2)
library(AER)
library(gridExtra)
library(estimatr)
library(summarytools)
library(tidyverse)
library(ggeffects)
library(sjPlot)
library(expss)
library(dotwhisker)
library(ggpubr)
library(broom)
library(dplyr)
library(plyr)
library(ggridges)
library(effects)
library(margins)
library(sandwich)
library(lavaan)
library(semTools)

load("Main Analyses/IndDatS4W3.RData")

mean(dat$apstot9, na.rm=T)
sd(dat$apstot9, na.rm=T)

mean(dat$otot9, na.rm=T)
sd(dat$otot9, na.rm=T)

mean(dat$atot9, na.rm=T)
sd(dat$atot9, na.rm=T)

mean(dat$mtot9, na.rm=T)
sd(dat$mtot9, na.rm=T)

## 9-item CFA
items<- dat[c("o1", "o6", "o8",
                "a1", "a2", "a4", 
                "m2", "m4", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o6 + o8
       a =~ a1 + a2 + a4
       m =~ m2 + m4 + m6
      '
#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", std.lv = T, items)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#N : 161
#robust CFI: 0.954
#robust RMSEA: 0.074
#latent cor: o/a-0.40, o/m-0.29, a/m- -.12
#omegas: O-0.69, A-0.80, M-0.849

#model setup 1 factor
cf1<- 'o =~ o1 + o6 + o8 +
        a1 + a2 + a4 +
        m2 + m4 + m6
      '
#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", std.lv = T, items)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI: 0.483
#robust RMSEA: 0.234

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

m1<- lm(ftrep ~ apstot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m1)

m2<- lm(ftdem ~ apstot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m2)

m3<- lm(ftparties ~ apstot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m3)

m4<- lm(idtot ~ apstot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m4)

m5<- lm(pvtot ~ apstot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m5)

m6<- lm(pktot ~ apstot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m6)

### subscales ###

m7<- lm(ftrep ~ otot9+atot9+mtot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m7)

m8<- lm(ftdem ~ otot9+atot9+mtot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m8)

m9<- lm(ftparties ~ otot9+atot9+mtot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m9)

m10<- lm(idtot ~ otot9+atot9+mtot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m10)

m11<- lm(pvtot ~ otot9+atot9+mtot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m11)

m12<- lm(pktot ~ otot9+atot9+mtot9+age+ba+white+hisplat+male+inc, 
       data = dat, na.action = na.exclude)

summary(m12)

gaze<- list(m1, m7, m2, m8, m3, m9)

stargazer(gaze,
          title = "",            
          dep.var.labels = "", 
          column.labels = c("ftrep", "ftdem", "ftparties"),
          column.separate = c(2,2,2),
          align = T,
          digits = 2,
          font.size = "scriptsize",
          style = "ajps",
          float.env = "table",
          keep.stat = "n",
          model.numbers = F,
          star.cutoffs = c(0.05, 0.01, .001))

gaze<- list(m4, m10, m5, m11, m6, m12)

stargazer(gaze,
          title = "Regression Results",            
          dep.var.labels = "", 
          column.labels = c("party social id", "gen pol violence", "pol knowledge"),
          column.separate = c(2,2,2),
          align = T,
          digits = 2,
          font.size = "scriptsize",
          style = "ajps",
          float.env = "table",
          keep.stat = "n",
          model.numbers = F,
          star.cutoffs = c(0.05, 0.01, .001))

